wineデータはワインの苦味を決める要因に関する実験データである。
考える要因は、ワイン生成時の気温(temp)の2水準(cold or warm)と
果肉と皮の接触の有無(contact)の2水準(no or yes)である。
9人の審査員(judge)が4通りある条件毎に2本ずつ(合計8本)の苦味を評価したデータである。
response : ワインの苦味を表すスコア(0-100点)rating : スコア(response)を5段階評価にグループ分けtemp : ワイン生成時の気温(cold or warm)contact : 果肉と皮の接触の有無(no or yes)bottle : ボトルの識別番号(合計8本)judge : 審査員の識別番号(合計9名)library(MASS)
data(wine, package = "ordinal")
str(wine)
## 'data.frame': 72 obs. of 6 variables:
## $ response: num 36 48 47 67 77 60 83 90 17 22 ...
## $ rating : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 2 3 3 4 4 4 5 5 1 2 ...
## $ temp : Factor w/ 2 levels "cold","warm": 1 1 1 1 2 2 2 2 1 1 ...
## $ contact : Factor w/ 2 levels "no","yes": 1 1 2 2 1 1 2 2 1 1 ...
## $ bottle : Factor w/ 8 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 1 2 ...
## $ judge : Factor w/ 9 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
\[ \log \frac{{\rm Pr}(rating \le j)}{{\rm Pr}(rating > j)} =\theta_j - \beta_1\cdot temp - \beta_2 \cdot contact \]
res <- polr(rating ~ temp + contact, data = wine, Hess = TRUE)
summary(res)
## Call:
## polr(formula = rating ~ temp + contact, data = wine, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## tempwarm 2.50 0.529 4.73
## contactyes 1.53 0.477 3.21
##
## Intercepts:
## Value Std. Error t value
## 1|2 -1.344 0.517 -2.600
## 2|3 1.251 0.438 2.857
## 3|4 3.467 0.598 5.800
## 4|5 5.006 0.731 6.850
##
## Residual Deviance: 172.98
## AIC: 184.98
polr関数の出力names(res)
## [1] "coefficients" "zeta" "deviance" "fitted.values"
## [5] "lev" "terms" "df.residual" "edf"
## [9] "n" "nobs" "call" "method"
## [13] "convergence" "niter" "lp" "Hessian"
## [17] "model" "contrasts" "xlevels"
回帰係数\( \beta_1,\beta_2 \)の推定値
res$coefficient
## tempwarm contactyes
## 2.503 1.528
カテゴリ間閾値\( \theta_1,\theta_2,\theta_3,\theta_4 \)の推定値
res$zeta
## 1|2 2|3 3|4 4|5
## -1.344 1.251 3.467 5.006
逸脱度 -2*logLik(res) と同じ
res$deviance
## [1] 173
カテゴリ選択確率 predict(res,type="probs") と同じ
res$fitted.value
## 1 2 3 4 5
## 1 0.206792 0.5706 0.1923 0.02362 0.006651
## 2 0.206792 0.5706 0.1923 0.02362 0.006651
## 3 0.053547 0.3776 0.4431 0.09582 0.029927
## 4 0.053547 0.3776 0.4431 0.09582 0.029927
## 5 0.020889 0.2014 0.5016 0.20049 0.075626
## 6 0.020889 0.2014 0.5016 0.20049 0.075626
## 7 0.004609 0.0538 0.3042 0.36359 0.273780
## 8 0.004609 0.0538 0.3042 0.36359 0.273780
## 9 0.206792 0.5706 0.1923 0.02362 0.006651
## 10 0.206792 0.5706 0.1923 0.02362 0.006651
## 11 0.053547 0.3776 0.4431 0.09582 0.029927
## 12 0.053547 0.3776 0.4431 0.09582 0.029927
## 13 0.020889 0.2014 0.5016 0.20049 0.075626
## 14 0.020889 0.2014 0.5016 0.20049 0.075626
## 15 0.004609 0.0538 0.3042 0.36359 0.273780
## 16 0.004609 0.0538 0.3042 0.36359 0.273780
## 17 0.206792 0.5706 0.1923 0.02362 0.006651
## 18 0.206792 0.5706 0.1923 0.02362 0.006651
## 19 0.053547 0.3776 0.4431 0.09582 0.029927
## 20 0.053547 0.3776 0.4431 0.09582 0.029927
## 21 0.020889 0.2014 0.5016 0.20049 0.075626
## 22 0.020889 0.2014 0.5016 0.20049 0.075626
## 23 0.004609 0.0538 0.3042 0.36359 0.273780
## 24 0.004609 0.0538 0.3042 0.36359 0.273780
## 25 0.206792 0.5706 0.1923 0.02362 0.006651
## 26 0.206792 0.5706 0.1923 0.02362 0.006651
## 27 0.053547 0.3776 0.4431 0.09582 0.029927
## 28 0.053547 0.3776 0.4431 0.09582 0.029927
## 29 0.020889 0.2014 0.5016 0.20049 0.075626
## 30 0.020889 0.2014 0.5016 0.20049 0.075626
## 31 0.004609 0.0538 0.3042 0.36359 0.273780
## 32 0.004609 0.0538 0.3042 0.36359 0.273780
## 33 0.206792 0.5706 0.1923 0.02362 0.006651
## 34 0.206792 0.5706 0.1923 0.02362 0.006651
## 35 0.053547 0.3776 0.4431 0.09582 0.029927
## 36 0.053547 0.3776 0.4431 0.09582 0.029927
## 37 0.020889 0.2014 0.5016 0.20049 0.075626
## 38 0.020889 0.2014 0.5016 0.20049 0.075626
## 39 0.004609 0.0538 0.3042 0.36359 0.273780
## 40 0.004609 0.0538 0.3042 0.36359 0.273780
## 41 0.206792 0.5706 0.1923 0.02362 0.006651
## 42 0.206792 0.5706 0.1923 0.02362 0.006651
## 43 0.053547 0.3776 0.4431 0.09582 0.029927
## 44 0.053547 0.3776 0.4431 0.09582 0.029927
## 45 0.020889 0.2014 0.5016 0.20049 0.075626
## 46 0.020889 0.2014 0.5016 0.20049 0.075626
## 47 0.004609 0.0538 0.3042 0.36359 0.273780
## 48 0.004609 0.0538 0.3042 0.36359 0.273780
## 49 0.206792 0.5706 0.1923 0.02362 0.006651
## 50 0.206792 0.5706 0.1923 0.02362 0.006651
## 51 0.053547 0.3776 0.4431 0.09582 0.029927
## 52 0.053547 0.3776 0.4431 0.09582 0.029927
## 53 0.020889 0.2014 0.5016 0.20049 0.075626
## 54 0.020889 0.2014 0.5016 0.20049 0.075626
## 55 0.004609 0.0538 0.3042 0.36359 0.273780
## 56 0.004609 0.0538 0.3042 0.36359 0.273780
## 57 0.206792 0.5706 0.1923 0.02362 0.006651
## 58 0.206792 0.5706 0.1923 0.02362 0.006651
## 59 0.053547 0.3776 0.4431 0.09582 0.029927
## 60 0.053547 0.3776 0.4431 0.09582 0.029927
## 61 0.020889 0.2014 0.5016 0.20049 0.075626
## 62 0.020889 0.2014 0.5016 0.20049 0.075626
## 63 0.004609 0.0538 0.3042 0.36359 0.273780
## 64 0.004609 0.0538 0.3042 0.36359 0.273780
## 65 0.206792 0.5706 0.1923 0.02362 0.006651
## 66 0.206792 0.5706 0.1923 0.02362 0.006651
## 67 0.053547 0.3776 0.4431 0.09582 0.029927
## 68 0.053547 0.3776 0.4431 0.09582 0.029927
## 69 0.020889 0.2014 0.5016 0.20049 0.075626
## 70 0.020889 0.2014 0.5016 0.20049 0.075626
## 71 0.004609 0.0538 0.3042 0.36359 0.273780
## 72 0.004609 0.0538 0.3042 0.36359 0.273780
線形予測子 model.matrix(res)[,-1] %*% res$coef と同じ
res$lp
## 1 2 3 4 5 6 7 8 9 10 11 12
## 0.000 0.000 1.528 1.528 2.503 2.503 4.031 4.031 0.000 0.000 1.528 1.528
## 13 14 15 16 17 18 19 20 21 22 23 24
## 2.503 2.503 4.031 4.031 0.000 0.000 1.528 1.528 2.503 2.503 4.031 4.031
## 25 26 27 28 29 30 31 32 33 34 35 36
## 0.000 0.000 1.528 1.528 2.503 2.503 4.031 4.031 0.000 0.000 1.528 1.528
## 37 38 39 40 41 42 43 44 45 46 47 48
## 2.503 2.503 4.031 4.031 0.000 0.000 1.528 1.528 2.503 2.503 4.031 4.031
## 49 50 51 52 53 54 55 56 57 58 59 60
## 0.000 0.000 1.528 1.528 2.503 2.503 4.031 4.031 0.000 0.000 1.528 1.528
## 61 62 63 64 65 66 67 68 69 70 71 72
## 2.503 2.503 4.031 4.031 0.000 0.000 1.528 1.528 2.503 2.503 4.031 4.031
\( (\beta_2,\beta_2,\theta_1,\ldots,\theta_4) \)推定量の分散共分散行列
vcov(res)
## tempwarm contactyes 1|2 2|3 3|4 4|5
## tempwarm 0.27950 0.04813 0.06543 0.12980 0.2240 0.2608
## contactyes 0.04813 0.22717 0.07483 0.12678 0.1730 0.2022
## 1|2 0.06543 0.07483 0.26739 0.09502 0.1022 0.1131
## 2|3 0.12980 0.12678 0.09502 0.19174 0.1936 0.2114
## 3|4 0.22403 0.17302 0.10218 0.19362 0.3573 0.3627
## 4|5 0.26084 0.20224 0.11315 0.21142 0.3627 0.5342
属するカテゴリの予測値
predict(res, type = "class")
## [1] 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3
## [36] 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3
## [71] 4 4
## Levels: 1 2 3 4 5
回帰係数の信頼区間
confint(res)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## tempwarm 1.5098 3.595
## contactyes 0.6158 2.492